#include <iostream>
#include <fstream>
#include <cmath>
#include <complex>
#include <cstdlib>
#include <vector>
using namespace std;

// 定义复数类型
typedef complex<double> Complex;

// 写入数据到文件
void writeDataToFile(const string& filename, const vector<Complex>& data) {
    ofstream file(filename);
    if (!file.is_open()) {
        cerr << "Error: Could not open file for writing." << endl;
        exit(1);
    }
    for (const auto& z : data) {
        file << z.real() << " " << z.imag() << endl;
    }
    file.close();
}

// Adams-Bashforth 方法
Complex bashforth1(Complex z) { return z - 1.0; }
Complex bashforth2(Complex z) { return z * (z - 1.0) / (1.5 * z - 0.5); }
Complex bashforth3(Complex z) { return z * z * (z - 1.0) / (23.0 / 12.0 * z * z - 16.0 / 12.0 * z + 5.0 / 12.0); }
Complex bashforth4(Complex z) { return z * z * z * (z - 1.0) / ((55.0 * pow(z, 3) - 59.0 * pow(z, 2) + 37.0 * z - 9.0) / 24.0); }

// Adams-Moulton 方法
Complex moulton3(Complex z) { return z * (z - 1.0) / (5.0 * z * z + 8.0 * z - 1.0) * 12.0; }
Complex moulton4(Complex z) { return z * z * (z - 1.0) / (9.0 * pow(z, 3) + 19.0 * pow(z, 2) - 5.0 * z + 1.0) * 24.0; }
Complex moulton5(Complex z) { return z * z * z * (z - 1.0) / (251.0 * pow(z, 4) + 646.0 * pow(z, 3) - 264.0 * pow(z, 2) + 106.0 * z - 19.0) * 720.0; }

// BDF 方法
Complex backward1(Complex z) { return (z - 1.0) / z; }
Complex backward2(Complex z) { return (3.0 * pow(z, 2) - 4.0 * z + 1.0) / (2.0 * pow(z, 2)); }
Complex backward3(Complex z) { return (11.0 * pow(z, 3) - 18.0 * pow(z, 2) + 9.0* z - 2.0) / (6.0 * pow(z, 3)); }
Complex backward4(Complex z) { return (25.0 * pow(z, 4) - 48.0 * pow(z, 3) + 36.0 * pow(z, 2) - 16.0 * z + 3.0) / (12.0 * pow(z, 4)); }

// 主函数
int main() {
    // 参数设置
    int numPoints = 2000;
    double thetaMax = 2 * M_PI;
    vector<Complex> zValues(numPoints);

    // 生成单位圆上的点
    for (int i = 0; i < numPoints; ++i) {
        double theta = i * thetaMax / numPoints;
        zValues[i] = polar(1.0, theta); // 单位圆上的点
    }

    // 计算并保存不同方法的稳定性函数值
    vector<Complex> bashforth1Data, bashforth2Data, bashforth3Data, bashforth4Data;
    vector<Complex> moulton3Data, moulton4Data, moulton5Data;
    vector<Complex> backward1Data, backward2Data, backward3Data, backward4Data;

    for (const auto& z : zValues) {
        bashforth1Data.push_back(bashforth1(z));
        bashforth2Data.push_back(bashforth2(z));
        bashforth3Data.push_back(bashforth3(z));
        bashforth4Data.push_back(bashforth4(z));

        moulton3Data.push_back(moulton3(z));
        moulton4Data.push_back(moulton4(z));
        moulton5Data.push_back(moulton5(z));

        backward1Data.push_back(backward1(z));
        backward2Data.push_back(backward2(z));
        backward3Data.push_back(backward3(z));
        backward4Data.push_back(backward4(z));
    }

    // 写入文件
    writeDataToFile("bashforth1.dat", bashforth1Data);
    writeDataToFile("bashforth2.dat", bashforth2Data);
    writeDataToFile("bashforth3.dat", bashforth3Data);
    writeDataToFile("bashforth4.dat", bashforth4Data);

    writeDataToFile("moulton3.dat", moulton3Data);
    writeDataToFile("moulton4.dat", moulton4Data);
    writeDataToFile("moulton5.dat", moulton5Data);

    writeDataToFile("backward1.dat", backward1Data);
    writeDataToFile("backward2.dat", backward2Data);
    writeDataToFile("backward3.dat", backward3Data);
    writeDataToFile("backward4.dat", backward4Data);

    // 使用 Gnuplot 绘图（Adams-Bashforth 修改部分）
    system("gnuplot -persist <<-EOF\n"
           "set title 'Stability Regions (Adams-Bashforth p=1, 2, 3)'\n"
           "set xlabel 'Re'\n"
           "set ylabel 'Im'\n"
           "set size square\n"
           "plot 'bashforth1.dat' with lines title 'Bashforth p=1', \\\n"
           "     'bashforth2.dat' with lines title 'Bashforth p=2', \\\n"
           "     'bashforth3.dat' with lines title 'Bashforth p=3'\n"
           "EOF");

    system("gnuplot -persist <<-EOF\n"
           "set title 'Stability Regions (Adams-Bashforth p=4)'\n"
           "set xlabel 'Re'\n"
           "set ylabel 'Im'\n"
           "set size square\n"
           "plot 'bashforth4.dat' with lines title 'Bashforth p=4'\n"
           "EOF");

    // Adams-Moulton 和 BDF 的绘图保持不变
    system("gnuplot -persist <<-EOF\n"
           "set title 'Stability Regions (Adams-Moulton)'\n"
           "set xlabel 'Re'\n"
           "set ylabel 'Im'\n"
           "set size square\n"
           "plot 'moulton3.dat' with lines title 'Moulton p=3', \\\n"
           "     'moulton4.dat' with lines title 'Moulton p=4', \\\n"
           "     'moulton5.dat' with lines title 'Moulton p=5'\n"
           "EOF");

    system("gnuplot -persist <<-EOF\n"
           "set title 'Stability Regions (BDF)'\n"
           "set xlabel 'Re'\n"
           "set ylabel 'Im'\n"
           "set size square\n"
           "plot 'backward1.dat' with lines title 'Backward p=1', \\\n"
           "     'backward2.dat' with lines title 'Backward p=2', \\\n"
           "     'backward3.dat' with lines title 'Backward p=3', \\\n"
           "     'backward4.dat' with lines title 'Backward p=4'\n"
           "EOF");

    cout << "Plots generated using Gnuplot. Check the output files and plots." << endl;

    return 0;
}